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Abstract 

A one-parameter family of second-order explicit and implicit total variation diminishing 
(TVD) schemes is reformulated so that a simplier and wider group of limiters is included. 
The resulting scheme can be viewed as a symmetrical algorithm with a variety of numerical 
dissipation terms that are designed for weak solutions of hyperbolic problems. This is a 
generalization of Roe and Davis’s recent works to a wider class of symmetric schemes other 
than Lax-Wendroff. The main properties of the present class of schemes are that they 
can be implicit, and when steady-state calculations are sought, the numerical solution is 
independent of the time step. 


I. Introduction 

The notion of total variation diminishing (TVD) schemes was introduced by Harten [1,2]. 
He derived a set of sufficient conditions which are very useful in checking or constructing 
second-order TVD schemes. The main mechanism that is currently in use for satisfying 
TVD sufficient conditions involves some kind of limiting procedure. There are generally 
two types of limiters: namely slope limiters [3] and flux limiters [4-6]. For a slope limiter 
one imposes constraints on the gradients of the dependent variables. In contrast, for a 
flux limiter one imposes constraints on the gradients of the flux functions. For constant 
coefficients, the two type of limiters are equivalent. The main property of a TVD scheme is 
that, unlike monotone schemes, it can be second-order accurate and is oscillation-free across 
discontinuities (when applied to nonlinear scalar hyperbolic conservation laws and constant 
coefficient hyperbolic systems). Sweby [5] and Roe [6] constructed a class of limiters as 
a function of the gradient ratio. Most of the current limiters in used are equivalent to 
members of this class. 

Although TVD schemes are designed for transient applications, they also have been ap- 
plied to steady-state problems [6-10], It is well known that explicit methods are usually 
easier to program and often require less storage than implicit methods, but can suffer a 
loss of efficiency when the time step is restricted by stability rather than accuracy. It is 
also commonly known that it is not very useful to extend Lax-Wendroff-type schemes to 
implicit methods, since the resulting schemes are not suitable for steady-state calculations. 
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This is due to the fact that the steady-state solution will depend on the time step. Roe 
has recently proposed a very enlightening generalized formulation of TVD Lax-Wendroff 
schemes [11], Roe’s results, in turn, is a generalization of Davis’s work [12]. It was the in- 
vestigation of these schemes which prompted the work of this paper. Their formulation has 
great potential for transient applications, but as it stands is not suitable for an extension 
to implicit methods. 

The aim of this paper is to incorporate the results of Roe [11], with minor modification, to 
a one-parameter family of explicit and implicit TVD schemes [2,8-9] so that a wider group 
of limiters can be represented in a general but rather simple form which is at the same time 
suitable for steady-state applications. The final scheme can be interpreted as a three-point, 
spatially central difference explicit or implicit scheme which has a whole variety of more 
rational numerical dissipation terms than the classical way of handling shock-capturing al- 
gorithms. In other words, it is a non-upwind TVD scheme or a symmetric TVD scheme. 
The proposed scheme can be used for time-accurate or steady-state calculations. Moreover, 
Roe’s formulation can be considered as a member of the explicit scheme for time-accurate 
calculations. By writing the scheme in terms of numerical fluxes with two input parameters 
(one for the choice of the time-differencing method and one for the option of choosing the 
Lax-Wendroff flux), a single computer program can be easily coded to include all of the 
schemes under discussion. Various limiters can be considered as external functions inside 
the computer program. Extension of the schemes to nonlinear scalar and system of hyper- 
bolic conservation laws is discussed in detail. An equivalent representation of the proposed 
numerical dissipation terms which avoids extra logic in the computer implementation is also 
included in the appendix. 

The effort involved in modifying some existing central difference computer codes for 
systems of hyperbolic conservation laws is fairly simple and straightforward. The proposed 
algorithms should have the potential of improving the robustness and accuracy of many 
practical physical and engineering applications. Only analytical formulations are presented 
here. Numerical testing and practical applications will be the subject of a forthcoming 
paper. 


II. Preliminaries 


In this section, a class of explicit and implicit TVD schemes [2] is reviewed. Harten’s 
sufficient conditions for this class of schemes are also stated. This set of conditions is then 
ultilized in the subsequent sections to construct and reformulate the second-order explicit 
and implicit TVD schemes of Harten [1,2]. 


Consider the scalar hyperbolic conservation law' 


du 

dt 


d/(tt) 
d x 


= 0 , 


( 2 . 1 ) 


where / is the flux and a{u) = df /du is the characteristic speed. Let u” be the numerical 
solution of (2.1) at x = jAx and t ~ nAt, with Ax the spatial mesh size and At the time 
step. Consider a one-parameter family of five-point difference schemes in conservation form 
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( 2 . 2 ) 


„“+> + - A”«) = «f - A(1 - 9)(/.” +i - 

where 0 is a nonnegative parameter, A = At/ Ax, ft n ±i — h(u^ 1 , u”, u” ±1 , u™.^), and 
«" +1 » ujnig). The function ft J+ i is commonly called a numerical 
flux function. Let 

fe J+ i=(l-^ + ,+^l (2.3) 

be another numerical flux function. Then (2.2) can be rewritten as 

u” +1 = u?-X{h j+ L -fty.i). (2.4) 

This numerical flux is a function of eight variables, ft - + jl = 

ft(w”_i, u™, u” +1 , u" +2 , u y-i lu > + 1 ’ u "+i> u y +2 )> anc * * s consistent with the conservation 
law(2.1) in the following sense 

h(u,u,u,u,u,u,u,u) = /(u). (2-5) 

This one-parameter family of schemes contains implicit as well as explicit schemes. When 
0 — 0, (2.2) is an explicit method. When 0^0, (2.2) is an implicit scheme. For example, if 
8 = 1/2, the time-differencing is the trapezoidal formula, and if 8 = 1, the time-differencing 
is the backward Euler method. To simplify the notation, rewrite equation (2.2) as 


L ■ u n+1 = R ■ u n , (2.6) 

where L and R are the following finite-difference operators: 

(L ■ u)y = Uj + A 0(ft 1+ j - hj_ j ) (2.7a) 

(#■«),- = Uj - A(l - 0)(ft y+ * - fty_i). (2.7b) 

The total variation of a mesh function u n is defined to be 

OO OO 

TV(u n ) - E |A i+4 «"|, (2.8) 

jzz — OO j— — OO 

where A J + xu n = u J + i — Uj. Here, the general notation convention 

Aj + iz — Zj-f-i — Zj (2-9) 

for any mesh function z is used. The numerical scheme (2.2) for an initial-value problem of 
(2.1) is said to be TVD if 

TV{u n +') <TV{u n ). (2.10) 
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The following sufficient conditions for (2.2) to be a TVD scheme are due to Harten [2]: 

TV (R ■ u n ) < TV (u n ) (2.11a) 

and 

TV{L-u n+l ) > TV{u n+1 ). (2.11b) 

Assume the numerical flux h in (2.2) is Lipschitz continuous and (2.2) can be written as 


n+l 

D 


u ; '--X6\C H ^A jJrk u- CJ iA^xu 


* + 

3-^3- 


n-f 1 


— u” 1 A( 1 0) , ; tl - c; I A.; l u j , 

( 2 . 12 ) 

where Cj^_ , = C+ (uj, u,'±i , ^± 2 ) or possibily C~ ±i = Uj,Uj±i,Uj± 2 ) ar e some 

bounded functions. Then Harten further showed that sufficient conditions for (2.11) are 

(a) if for all j 


C? 4J =A(1-*)C£ 4 >0 


c+ +! +c; 4i = mi - »)(c; +J +c; +4 ) < 1, 


(2.13a) 

(2.13b) 


and 

(b) if for all j 


< C < XOC 




3+i 


< 0 


(2.14) 


for some finite C. Conditions (2.13) and (2.14) are very useful in guiding the construction of 
second-order accurate TVD schemes which do not exhibit the spurious oscillation associated 
with the more classical second-order schemes. 

Harten [1-2], Yee et al. and Yee 7-10] investigated a particular form of C ± . They 
have shown in a variety of numerical tests that the scheme is quite useful for gas-dynamic 
calculations. The recent work of Roe [ill suggests a wider class of flux limiters for the 
Lax-Wendroff-type of TVD schemes which with a minor modification are found to have 
an immediate application to scheme (2.2). The details will be discussed in the next two 
sections. 


III. A Generalized Formulation of a Class of Symmetric Schemes 

In this section, Roe’s formulation is reviewed. Then, with a minor modification, his 
numerical flux is shown to be applicable to a larger class of symmetric schemes. Sufficient 
conditions for this new class of schemes to be TVD are derived for both the constant 
coefficient and nonlinear scalar hyperbolic equations. 
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3.1 Roe’s TVD Lax-Wendroff Schemes 


Roe [11] has recently developed a generalized formulation of TVD Lax-Wendroff schemes. 
The form of the schemes is the usual Lax-Wendroff plus a general conservative dissipation 
term designed in such a way that the final scheme is TVD. For 5/ jdu = a — constant, his 
scheme is written as 


M y +1 = U 1 _ \ v ( l + v )^3-h u ~ ~ 

“ ^H( x - MK 1 

+ ^M(l- W\){1-Q 3+ i)A^xu. (3.1) 

Here v — aX = aAt/ Ax, the first two terms represent the usual Lax-Wendroff scheme, and 
the other two terms represent an additional conservative dissipation. The function Qj+i 
depends on three consecutive gradients A A j + iu, A J+ aii and is of the form 

Q j+ ±=Q(rJ+ i,r+ + i), (3.2a) 

where 

A y_JLtt Ay + |« 

Ay +J LU’ ^'+2 Ay + iU 

Here are not defined if A + j u — 0. To avoid this, an equivalent representation will 
be discussed in the appendix. If one assumes both Q and Q/r are always positive, then a 
set of sufficient conditions for (3.1) to be TVD is 


(3.2b) 



Qj+i < 

2 

(3.3a) 


1 - M 

Qj+I j 


2 

|j/| 

(3.3b) 

Qj+ 5 / 

/-m)< 

2 

H' 

(3.3c) 


Two examples for the function Q are 

Q(r~ ,r + ) = minmod(l, r~) + minmod(l, r + ) - 1, (3-4) 

and 

Q(r~ ,r*) = minmod(l,r'”,r + ). (35) 

Normally the “minmod” function of two arguments is defined as 
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minmod(a:, y) = sgn(x) • max{0,min[|x|,y • sgn(x)] }. 
but within this context 

j/i i'* Jmin(l,r ± ) r* > 0 /0 

minmod^r*) = | Q v ' r ± <- Q (3-6) 

Other forms of Q(r ~ , r + ) are discussed in Sweby [5] and Roe [6], 

Scheme (3.1) is a reformulation of Davis’s work [12] in a way which is easier to analyze 
and includes a class of TVD schemes not observed by Davis. The numerical flux denoted 
by h h ™ L for (3.1) is 

^•■,4 r = 2 { a ( w i+i + u j) ~ [^ a 5 + MU — • (3-7) 

Scheme (3.1) is second-order accurate in time and space. Observe that by setting 9 — 0 in 
(2.2) and by using (3.7) as the numerical flux, the resulting scheme is (3.1). 

3.2 Schemes for Linear Scalar Hyperbolic Equations 

If one is to use (3.7) as the numerical flux for (2.2) with 9^0, then the resulting scheme 
is only useful for transient calculations. For steady-state applications, either one has to 
restrict the time step similar to the explicit method or the steady-state solution will depend 
on the time step. It is emphasized here that the dependence of the time step in steady-state 
solutions occurs even though the value of At is similar to an explicit method. In this case 
At is most often in the same order as Ax; thus the dependence in At is less severe. The 
term that causes this undesirable property is the one with coefficient a 2 in equation (3.7), 
since the A appears in this term. Therefore, besides considering the use of (3.7) as the 
numerical flux for (2.2) when 9 = 0, the numerical flux (3.7) with a 2 — 0 is also considered; 
i.e.. the numerical flux is of the form 

«(«y+ 1 + u i) ~ NU - ■ (3.8) 

Now the question is, will the new numerical flux (3.8) satisfy the sufficient conditions (2.11)? 
The answer is yes. It turns out that some of the Q functions that are suitable for the 
generalized TVD Lax-Wendroff scheme are also suitable for (3.8). The implication is that 
if one chose the proper Q function, the resulting scheme (2.2) together with (3.8) can be 
viewed as a symmetrical algorithm with a wide variety of numerical dissipation terms that 
satisfy the TVD property. 

Now with the choice of (3.8), the corresponding C ± of equation (2.12) are 

~ 2^-2 2 + ’ a > 0 (3.9a) 



6 



a < 0 . 


(3.9b) 


C ',H = 101 


1 2® , + i + 2 r >~i 


Therefore, sufficient conditions for this specific numerical flux function (3.8) to be TVD are 

1 “ + 2 (^ H 


0 < A (1 — 9)a 
0 < A (1 - 0)\a\ 


and 


Q >+i/ r i*i) 

< 1 

a > 0 

(3.10a) 


< 1 

a < 0 

(3.10b) 

-A0C± < 0. 

J + 5 ~ 



(3.11) 


For 0 < 6 < 1 and v 7 ^ 0, condition (3.10a) is satisfied if 


^■+i/ r 7+4 

and condition (3.10b) is satisfied if 


't-H I <2 


Qj-i < 
A a < 


A(l - 9)a 

1 


- 2 


1 - e ’ 


[Qi-i/rt, <2, 




J-2 


i / < 


J + 2 A(1 - 0)|a| 

1 


Ala < 


\-e ' 


(3.12a) 

(3.12b) 

(3.12c) 


(3.12d) 
(3.12e) 
(3.1 2f) 


Since A|a| < — j, the term - 2 is always positive. Therefore the same assumption 

as Roe can be made; i.e., assume both Q and Q/r are always positive. Then all one has to 
do is devise a function Q such that 


Q 0+ i < 2 (3.13a) 

(o*ti/-7+*) < 3t(r^H -* (313b) 

(«. + i/ r ;o)<X(T^jR - 2 ( 313c ) 
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.1 I 1 

A o < 

1 1 1-9 


(3.13d) 


With the above choice of Q, the last sufficient condition (3.11) is immediately satisfied. 
Some readily available Qj + i functions can be found in references [5,11]. For instance, the 
two examples given in equations (3.4) and (3.5) satisfy condition (3.13). 

For 9 — 1/2, scheme (2.2) together with (3.8) and (3.13) is second-order accurate in 
both space and time. The CFL-like restriction for (2.2) to be TVD in this case is 2. When 
9 — 1, scheme (2.2) together with (3.8) and (3.13) is unconditionally TVD, but the resulting 
scheme is first-order in time and second-order in space. When 0 = 0, the scheme is explicit, 
and unlike Roe’s schemes, is only first-order in time but second-order in space. 

As noted before, the value , , (or r + , ,) is not defined if A. iu (or A =tf) is finite 
but Aj + iu = 0. For computer implementation purposes, it might be more convenient to 

define Qj + i A J+ iu = + where Q J + ± is a function of A y_i«, Ay + iu, and A ?+ 3 u, but 

not a ratio of those gradients. For this formulation and its extension to the nonlinear case, 
see the appendix. 

3.3 Linearized Version of the Proposed Scheme 
for Constant Coefficient Equations 

For 0 -■£ 0, scheme (2.2) is implicit. Moreover, this is a genuinely nonlinear scheme in 
the sense that the final algorithm is nonlinear even for the constant coefficient case. The 
value of u n+1 is obtained as the solution of a system of nonlinear algebraic equations. To 
solve this set of nonlinear equations noniteratively, a linearized version of (2.2) together 
with (3.8) is considered. Substituting (3.8) in (2.2), one obtains 


u] +1 + ~ 
3 2 

A 0 

~ Y 


auj+i ju|(l ■ i ) A j_|_ i u 

i - M(1 - <9y-i)A 


n+ 1 


au 


lU 


n+ 1 


= RHS of (2.2). (3.14) 


Here “RHS of (2.2)” means the right hand side of equation (2.2) with h J + 1 defined in (3.8). 

Locally linearizing the coefficients of (A J± iu) n+1 in (3.14) by dropping the time index from 
[n + 1) to n, one gets 


n + 1 


A 9 

Y 


au 


n-f 1 

i+ 1 ' 


N(i - e" + j)A, + j 


U 


n-f 1 


o|(l — Q n _ i )A •_ i it n+1 = RHS of (2.2). (3.15) 


Lst d j 


u” +1 — u”; i.e., the “delta” notation, equation (3.15) can be written as 
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(3.16) 


eidj-i + ezdj + e 3 < f J+1 — — A(A" + | — ft"_±) 

where 


A 0 


ei 


a ~ |a|(l-Q,_l) 


e 2 = 1 + 


xe 


|o|(i-g i) + |o|(i-g, + i) 


xe 


e 3 


a - |a|(l-g jH . 


(3.16b) 

(3.16c) 

(3.16d) 


The linearized form (3.16) is a spatially five-point scheme and yet it is a tridiagonal system 
of linear equations. This is because at the (n+l)th time level, only three points are involved; 
i.e., u^l,u^ + \u^l Although the coefficients e % involve five points, they are at the nth 
time level. 

The form of , for (3.15) is the same as (3.9) except the time index for the Qj±± and 
, is dropped from (n-H 1) to n for the implicit operator. Thus the linearized form (3.16) 

is still TVD. It was found in reference [7] that when time-accurate TVD schemes are used 
as a relaxation method for steady-state calculations, the convergence rate is degraded if 
limiters are present on the implicit operator. For steady-state applications, one can obtain 
another TVD linearized form by setting Qj±i — 0 in (3.16); i.e., by redefining (3.16) by 


M, , o 

e i = y (- fl “ l fl l) 

(3.17a) 

e 2 — 1 + A0 (|a|) 

(3.17b) 

i 

il 

CO 

QA 

(3.17c) 


Scheme (3.16a) together with (3.17) is spatially first-order accurate for the implicit operator 
and spatially second-order accurate for the explicit operator. It can be shown that (3.16a) 
together with (3.17) is still TVD. Equation (3.17) is considered because no limiter is present 
for the implicit operator. 


3.4 Scheme for Nonlinear Scalar Hyperbolic Conservation Laws 


To extend the scheme to nonlinear scalar problems, one simply defines a local character- 
istic speed 


V-H 


( d f/ du ) L 


A y + iU ^ 0 

a m« = ° 


(3.18) 
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(3.19) 


and redefines the r^ +A in (3.2b) as 


r j+i 


K-xlAj.-ru 

K+§l A j+* tt ’ 


r + 

3+2 


l a j+|l A 3+f U 


l°i+*l A i+f tt 


Unlike the constant coefficient case, dj+i and a 3 -_ i are not always of the same sign. After 
considering all the possible combinations of the signs of the a J+ 1 and Oy_ i , a set of sufficient 
conditions on Q still can be of similar form to (3.13) and is 


Qj+i < ^ 


Q, 


.■+i/ r ,-+* 1 < 


A(1 - 0)\a-_; 




r 


< 


j+ *J A(l-fl)|o i+ a 

A|0 > + si < (T^)' 


- 2 


(3.20a) 

(3.20b) 

(3.20c) 

(3.20d) 


The numerical flux for the nonlinear case is 


= 2 


if 3 + 1 + Ji) ~ K+f K 1 ~ 


(3.21) 


Observe that when a - + ± = 0, the scheme has zero dissipation. One way is to approximate 
|n J+ i| by a Lipschitz continuous function [2], For example, instead of using (3.21), one can 
use 


h >*i = \ 


if 3+1 + fj) - X 1 - Qi+0 A s+i u 


Here ip is & function of a J+ x and is of the form 


or 


tp(z) = 


1 p{z) 


\z\ > e 


(z 2 + t 2 )/ 2e |z| < ( 


| z\ > e 
z < e 


where c is a positive small number [9], 


(3.22a) 


(3.23) 


(3.24) 


3.5 Alternate Scheme for the Nonlinear Scalar Hyperbolic Problem 

A simpler but less rigorous way of extending the constant coefficient case to the nonlinear 
case is to define a local characteristic speed j and keep the restriction on Q the same as 
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in (3.18) and (3.20), but use the . in (3.2b) instead of (3.19). The alternate form requires 

* J ' 2 

less computation than the previous more rigorous approach. The relative advantage and 
disadvantage between these two forms remain to be shown. However, numerical experiments 
with two-dimensional Euler equations of gas dynamics [7,9] show that the alternate form 
gives a better shock resolution than the former one (3.18)-(3.22). 

As a side remark, a case of Harten’s second-order explicit TVD scheme is contained in 
the class of limiters of Sweby [5] and Roe [6] and is equivalent to a case of Roe’s second- 
order scheme of reference [6], See Sweby ’s original manuscript [13] instead of the published 
version [5] for details. The numerical experiments of Yee et al. and Yee [7-9] with Harten’s 
second-order TVD scheme indicate that the alternate form is favored over the more rig- 
orous approach (3.18)-(3.22). This indication is further endorsed by Davis’s numerical 
experiments [12] with similar examples. Since different limiters have different effects which 
are highly problem-dependent on the resolution of the numerical solutions, all these possi- 
bilities in extending to nonlinear equations require more extensive numerical testing before 
a clearer picture can be drawn. 


3.6 Linearized Version of the Proposed Implicit Scheme for Nonlinear Equations 

For the nonlinear case, the situation is slightly more complicated since the characteristic 
speed df jdu is no longer a constant. Substituting (3.22) in (2.2), one obtains 


«r 1 + t 


\6 

2 


fj+ 1 V , ( a j+ 

fj - 1 - V’K- 


){l Qj+h)^J+% U 
)(1 - 


n+1 


n + 1 


= RHS of (2.2). 


(3.25) 


Here “RHS of (2.2)” means the right hand side of equation (2.2) with h } + \ defined in (3.22). 
Unlike the constant coefficient case, one also has to linearized j 1 , and Q"±l- 

Following the same procedure as in [9], two linearized versions of (3.25) are considered. 


Linearized Nonconservative Implicit Form 

Add and substract /” + 1 on the left-hand-side of (3.25) and using the relation (3.18), one 
can rewrite (3.25) as 


u“ +1 + - 
^ 2 

X0 

2 


a n+ i-^i)(i-Q: +u 




A+U 


1 u 


n+ 1 




A J _iu n+1 = RHS of (2.2). (3.26) 


By dropping the time index of the coefficients of Ay±iu n+1 from (n + 1) to n and letting 
u", (3.26) becomes 


di = < +1 
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(3.27a) 


eidj-i + e 2 dj + ezdj+i — -A {h* + i — A"_a) 

where 


and 


ej = A 6B 

e 2 = 1 - \o(b~ + £ + 
e 3 = X9B+ 


(3.27b) 

(3.27c) 

(3.27d) 


B~ 


±a j±k - ^(a i± i)(l - Qj±i) 


(3.27e) 


Again equation (3.27) is a five-point scheme, and yet the coefficient matrix associated 
with the dj ’ s is tridiagonal. With this linearization, the method is no longer conservative. 

^ - 4 - 

Therefore (3.27) is only applicable for steady-state calculations. The form of C ~_ ( _ 1 (which 

is the counterpart of (3.9) for the nonlinear case) is undisturbed, except the time index 
is dropped from (n + l) to n for the implicit operator. Thus the linearized form (3.27) 
is still TVD. Again, a spatially first-order accurate implicit operator similar to (3.17) can 
be obtained for (3.27) by setting B* = ^[±a J± i - ^(ia^ a)]”. Since the limiter does 
not appear on the left-hand-side, improvement in efficiency over (3.17) might be possible 
[7,9], This reduced form is especially useful for multidimensional, nonlinear, hyperbolic 
conservation laws. 


Linearized Conservative Implicit Form 

One can obtain a linearized conservative implicit form by using a local Taylor expansion 
about u n and expressing f n+1 - f n in the following form 

/; +1 - /; = a i ( u y 4 1 “ O + 0(At 2 ), (3.28) 

where a" = ( df /du )”. Applying the first-order approximation of (3.28) and locally lin- 
earizing the coefficients of (A J± iit) n+1 in (3.25) by dropping the time index from (n + 1) 
to n, one gets 


u n + 1 + — 
3 2 


“3-^-1 - V’K+rXl - Q ” + x)A m 


n-f 1 


+V’K_J(i-Q"_ i )A, 1 


n+1 


RHS of (2.2). (3.29) 


Letting d 3 = u 


n + l 


a", equation (3.29) can be written as 
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e i dj-i + e%dj + e^dj+i — — A(/i"_|_± — h"_i)> (3.30a) 

where 

(3.30b) 
(3.30c) 
(3.30d) 


xe 


ei 




e 2 = 1 + 


A0 


V’K-iM i - Qj_i) + V’K+i'K 1 - Qj 


j+i 


xe 


e 3 = 


«j+i - V'(s+i)(i - Qj+i) 


The linearized form (3.30) is conservative and is a spatially five-point scheme with a tridiag- 
onal system of linear equations. Scheme (3.30) is applicable to transient as well as steady- 
state calculation. But the form of for (3.30) is no longer the same as its nonlinear 

counter part. As of this writing, the conservative linearized form (3.30) has not been proven 
to be TVD. 

For steady-state application, one can use a spatially first-order implicit operator for (3.30) 
by simply setting all the Qj±i = 0; i.e., redefine (3.30b)-(3.30d) as 


X9 


e\ 


e 2 = 1 + 


xe 


es 


-Oj-l - 

xe [ 

+ V'K'+j) 


a t+i 


(3.31a) 

(3.31b) 

(3.31c) 


In reference [7,9], this type of linearization proved to be very useful for two-dimensional 
steady-state airfoil calculations. 


IV. Extension to Hyperbolic System of Conservation Laws 

Extension of the scalar scheme (3.14), (3.27), or (3.30) to systems of conservation laws 
can be accomplished by defining at each point a “local” system of characteristic fields, and 
then applying the scheme to each of the m scalar characteristic equations. Here m is the 
dimension of the hyperbolic system. Extension of the scalar scheme to higher than one- 
dimensional systems of conservation laws (for practical calculations) can be accomplished 
by an alternating direction implicit (ADI) method similar to the one described in Yee et al. 
and Yee [7,9], Only the one-dimensional case will be described here. 


Formal Extension 
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Consider a system of hyperbolic conservation laws 


w W| 

dt dx 


(4.1) 


Here U and F(U) are column vectors of m components. Let A = dF/dU and the eigenvalues 
of A be (a 1 , a 2 , Denote R {R~ 1 ) as the matrices whose columns are right (left) 

eigenvectors of A (yW 1 ). Let (7y + i denote some symmetric average of Uj and U J+1 (see 

references [1,7,9] for a formula). Let al J + i, + R~+i denote the quantities a 1 , R, R ~ 1 


evaluated at f/ J + i . Define 


j + i 


"j + 4 = R 3 + i A J + i U 


(4.2) 


as the forward difference (or the jump) of the local characteristic variables. With the above 
notation, a one-parameter family of TVD schemes (2.2) in the system case can be written 
as 


U 


n-f 1 


+ - h^D = \ - A(i - e)(H* +i -#»_,). 


rn - f 1 


i-H 


The numerical flux function H J+ 1 expressed as 


H i+i = 2 + F ’ +1 " J ’ 

where the elements of the denoted by <f> l j + ±,l — 1 are 

= V'(a'+i)(i-gi + iK + i, 

with t/)(^) defined in (3.23), and 


(r; + ,)^(r; +i )* 






\a[ ! | a' , 

1 J- 5 1 J- 2 


a ■ 


a 


j+ 2 y+5 


iV§l a j+§ 

K + 4 K +4 


(4.3a) 

(4.3b) 

(4.3c) 

(4.3d) 

(4.4a) 

(4.4b) 


or 


<v 4 >' 






Q: 


f+2 




(4.4c) 

(4.4d) 
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Here ct l . ± are the elements of (4.2). The corresponding conservative linearized form (3.30) 
2 

for the system case can be expressed as 

EiDj - 1 + EiDj + E 3 D H1 = -A (tf n +| - tf”_±), (4.5a) 

where 


with 


e ' = T 

E2 = I+ y(K,_s + K J + 4 ) 

\o f \ n 

E s = T {^A j+ 1 -K J+ ,J , 


and 



R-'QR 



n 


j ± ; 


diag 


^(a' ±i )(l-Q' ±i ) 


J±1 


or for the first-order left-hand-side 


n j± j = diag 


V’Kir) 


(4.5b) 


(4.5c) 

(4.5d) 


(4.5e) 


(4.5f) 


(4.5g) 


Here diag(z { ) denotes a diagonal matrix with diagonal elements z l . Aside from computing 
the right hand side, the rest of the arithmetic involved for equation (4.5) is two matrix mul- 
tiplications and a block tridiagonal inversion. The value of O + i in (4.5f) or (4.5g) can be 
saved while calculating the right hand side. Similarly, one can express the nonconservative 
linearized form (3.27) for the system case. 

As a side remark, with the same procedure Roe’s numerical flux in system case can be 
written 


HW = 






: 3 + k 


(4.6a) 


/ sLW 

where the elements of the denoted by ( < f >l ] + i) J — 1, are 


_ + a j+i 


(4.6b) 
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Simplified Version 

As one can see, the main work for scheme (4.3a) with numerical flux function (4.3b) 
or (4.6) is the term + A similar situation is also true for the corresponding 

conservative or nonconservative linearized form. Since 


H 7+i A <+ i v < ( 4 - 7 > 

if somehow one can simplify RUR~ X to be a diagonal matrix, then the current implicit 
scheme should be competitive in terms of operation count with the widely distributed codes 
such as ARC2D (version 150) of Pulliam and Steger [14] and FL052R of Jameson et al. 
[15], Both ARC2D and FL052R use a spatially three-point central differencing scheme 
with identical numerical dissipation terms, but they use different time-stepping methods 
for steady-state calculations. 

Two possible ways of simplifying (4.7) are one suggested by Davis [12], and one suggested 
by Roe [11]. Davis suggests approximating (4.7) with 


R 3+h®i+h = fi j+i n i+4 a j+i 


R i+^i- 


iU" 1 , A 

2 3+2 


3+i U 




diag 


V'K+i )(i 


Q l 


3 + h 


R ■ , id.'ii R 


l _ — 


v+rD+r j+i - ~^j+r'3+ 




(4.8) 


where — a; ( r ,^i5 r ^_ i ) ls a scalar function of an ^ a \ + ±-‘ The s y m bol I i 


3 ' H 




3 + ~r 


3+$ 


y+ i 


IS 


a rn X m identity matrix. Adapted to the current scheme, c can be expressed as 


“j+i = ^( a 3+0 


1 - Q{r : + i,r+ +i ) 


j+i 


(4.9) 


where 


a i+h 


max Jl a J + i |. 


(4.10) 


In order to not have to compute R and R 1 , have to be redefined such that they are 
functions of gradients of the original variables U instead of a • , i. Davis suggests using 


(4.11a) 




Ei(^i-i«')( Ay+i***) 

r 3+-k 

(A,. + .£7,A jH U) - 

Ei(A J+ iu') 2 

r + 

(A J+i t/,A,. + jt7) 

Ei(A i+ |«')( A y+i u ') 

?+! 

(A J+| f/,A iH U) - 

£«( a,- +4 «0 8 


(4.11b) 


where (•,•) denotes the usual inner product on the real m-dimensional vector space R m . 

Roe suggests that instead of using the fastest wave (4.10), one might consider the use of 
the strongest wave in the following sense 
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(4.12a) 


- _ zr=y (« 1 ) 2 

E i = i (« i ) 2 


i+i 


or 


- EZii^y 

J+ * zr=i(° i ) 2 


j+i 


(4.12b) 


For the one-dimensional Euler equation of gas dynamics (perfect gas), the a 1 are simply 


a 1 = u - c, a 2 = u, a 3 = u + c, 


(4.13) 


and a formula for the a 1 are 


„ i 

a - , i 


* 5 +1 


A j+i p ~ pj+j c j+jAj+V 


„ 2 

°J+i 


r 2 \ C i+i A i+5^ 

c i+4 


a 


,+ * * 5+4 


— A ,'+if ,+ '\+i c J+i A J+i“ ■ 


(4.14a) 

(4.14b) 

(4.14c) 


Here u is the velocity, c is the sound speed, p is the pressure, and p is the density. 

As for the r ± , Roe suggests the use of a J± ± instead of A j±iU in (4.11). In this case, 
the operation count between the local characteristic variable approach (4.5) and Roe’s 
suggestion might be very competitive, since in Roe’s suggestion one has to compute all 
the a ( . + A anyway. Therefore, the main difference is between computing = 


R J + idiag ij){a l . ±L ){ 1 - Q l j±i) together with (4.4), or computing equations (4.8)-(4.9) and 
equation (4.12a) or (4.12b) together with 


(ttj-i.a j+ i) 


(4.15a) 

r ;’+ 2 

(a i+ £,a i+ |) 

E,K +I ) 2 

r + 


_ ^' Q i+§ a ‘i+$ 

(4.15b) 

j+i 

K+b a y+f) 

E,(«J +i ) 2 ' 


In summary, for the system case three approaches are suggested here, (a) The more 
systematic approach (4.5) (from here on referred to as the local characteristic approach), 
(b) Davis’s approach, and (c) Roe’s approach. Davis’s suggestion is by far the simpliest 
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to implement and requires the least operation count but is the least rigorous. In numeri- 
cal experiments with a two-dimensional shock reflection problem and a circular arc airfoil 
problem, Davis’s approach shows good potential. Roe’s suggestion, in the author’s opinion, 
seems more rigorious than Davis’s, but requires a similar operation count to the current 
suggestion. An obvious advantage of the Davis and Roe approaches over the local charac- 
teristic approach is that one does not have to deal with the problem of singular values of 
and in turn Qj±i are rarely zero. An analytical, but not necessarily practical, advan- 
tage of the current aproach over the other two approaches is that (4.5) collaspes into the 
exact scalar scheme for the case m = 1. The implication is that if one locally freezes the 
coefficients in (4.5), then the resulting constant coefficient system is TVD and convergent 
subject to the CFL restriction of (frjj- The proof of this statement is readily available in 
reference [2], The total variation definition for a vector grid function for the system case is 

m 

Tv<m = X2>; + ji («e) 

3 i = i 


V. Concluding Remarks 

The present paper was inspired by the work of Roe [11] and Davis [12], and is based on the 
work ofHarten [1-2] and of Harten and the author [7-10]. A one-parameter family of explicit 
and implicit TVD schemes is reformulated so that a wider group of limiters is included. The 
current class of schemes as well as Roe and Davis’s can be classified as non-upwind TVD 
schemes or symmetric TVD schemes. The main advantages of the present class of schemes 
over the ones suggested by Osher and Chakravarthy [16], Roe, or Davis are that (a) a wider 
class of time-differencing is included, (b) the implicit scheme allows a natural linearized 
procedure for a noniterative implicit procedure, and thus might have a greater potential for 
practical applications, especially for “.stiff” problems, and (c) when applied to steady-state 
calculations, the numerical solution is independent of the time step. Furthermore, Roe’s 
formulation can be considered as a member of this family by simply setting 0 = 0 and using 
the numerical fluxes (3.7). Extension of this class of schemes and Roe’s schemes to a system 
of equations is straightforward. One can define a general numerical flux function 




- Rj , 1 < 1 > , . 


(5.1a) 


where the elements of the 4*^ , denoted by (</F + 4 ) ,/ = 1 , ...,m 

' = W{a l j+i) Qj+k + K+iK 1 ~ 


are 


tt i+ 4 ’ 


(5.1b) 


with P — 1 or 0 for transient calculations, and ft = 0 for steady-state calculations. Here 
when p — 0, (5.1) is (4.3b), and when p — 1 , (5.1) is the Roe’s Lax-Wendroff numerical 
flux (4.6). The suggestions of Roe and Davis to simplify the system case are also discussed. 
Numerical testing and practical applications will be carried out in a separate paper. 
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The results of Roe, Davis, and the present formulation provide a more rational way of 
supplying additional numerical dissipation terms to the commonly known schemes such 
as the Lax-Wendroff type and some spatially symmetrical explicit and implicit types of 
schemes. Here the amount of work required to modify existing computer codes with the 
suggested numerical dissipation terms varies from very minor changes to moderate yet 
straightforward computer programming. The potential of improving the robustness and 
accuracy of a wide variety of physical applications is worth the effort of further pursuing 
the implementation of these ideas into the many existing user-oriented computer codes. 


Appendix 

(Equivalent Representation for the Conservative Dissipation Term) 

The terms r ± t , in (3.2b) are not defined if A,_i« and A, + j u are finite and A.iiu — 0. 

3-r 2 3 2 JT 2 J ~ 2 

To avoid the use of extra logic in a computer implementation, it might be better to rewrite 
the terms Q ]+ i iu in equations (3.1), (3.8), and thereafter in the form 




(A.l) 


Linear Scalar Hyperbolic Equations 

The form Qj + i is a function of Ay_j.it, A J + j u, and Ay + att, but not r ±, ) i.e., 

Qj+i ~ Q(Ay_iU,A y+ iu,A 3 + |it). 

The numerical flux /iy + i in (3.8) can be rewritten as 

1 [ ^ 

Q ( Uj y j -f tty) — |n | ( Aj_|_ i It — Q y_y i 


IV • I 1 

r+2 2 

The corresponding C ± in (3.9) becomes 


C + x = a 1 - 


l Q , 


+ 


l Q 


j+i 


2 A j-iu 2Aj_iu/ : 


l <5, 


1 Qj--; 
+ -- ■ 


C _ , , = la 1 - 

3+2 V 2 1 U 2 A J+ iU 

The sufficient conditions for TVD in terms of Q are 


a > 0 


a < 0. 


(A.2) 


(A.3) 


(A. 4a) 
(A. 4b) 


A,-_ i u 


Q 


2 


A 3-± U 


< 2 


(A. 5a) 
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for a > 0, and 


gjf+j _ Q±± < 2 

Ay_iU Ay_iu A(l - 0)a 

1 


a < 


1-6 




A J+i u 


Qj- 


< 2 




< 


A J+ i u A J+ i« A(1 - 0)|a| 


a < 


1 - 


(A. 5b) 
(A. 5c) 


(A. 6a) 

(A. 6b) 
(A. 6c) 


for a < 0. Assume j + x / &j + xuj and are always positive. Then 

sufficient conditions for (A. 5) and (A.6) are 


Qj±i_ 

A H± U 

Q 


< 2 , 


]+ i 


< 


A •_ xu A(1 — 0)|a| 


J 2 
''3 


Q i— i 2 


< 


A j + iu A(1 - 0)|a| 


- 2 . 


(A. 7a) 
(A.7b) 
(A.7c) 


The expressions (3.4)-(3.5) now become 


^(Aj-i u > Ay+i“,A y+ au) = minmod (a j+ iu, A y _ii 


+ minmod I Ay + 1 u, A 3 + |U | -A J+ iu, 


(A. 8a) 


and 


Q [ A 7 -_iU,A i+ JLU,Ay + iU 


= minmod ^Ay_±u, Ay + iu, Ay + a« 


(A. 8b 


Here 
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(A.9) 


minmod(z, y) = sgn(:r) • max{0,min [|a:|, y • sgn(z)] }. 

In general, the “minmod” function of a list of arguments is equal to the smallest number in 
absolute value if the list of arguments is of the same sign, or is equal to zero if any argument 
is negative. 

Nonlinear Scalar Hyperbolic Conservation Laws 

For nonlinear problems, one way is to replace all the a’s in equation (A.3)-(A,7) by a J± i 

accordingly. The value of a - + 1 is defined in (3.14). For a more rigorous approach Q should 
also be redefined as 

Qj+$ = Q(l«i-^|A i _AU,|ay + i| A J - + iu,| a J - + !|A i+ au), (A.10) 

and A j±±u should be replaced by \aj±i |A J± i u whenever they appear in equations (A.4)- 
(A.9). Similarly, the system case can be rewritten in terms of the Q’s. 
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